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Abstract. We consider the exact path sampling of the squared Bessel process and some other 
continuous-time Markov processes, such as the CIR model, constant elasticity of variance dif- 
fusion model, and hypergeometric diffusions, which can all be obtained from a squared Bessel 
process by using a change of variable, time and scale transformation, and/or change of mea- 
sure. All these diffusions are broadly used in mathematical finance for modelling asset prices, 
market indices, and interest rates. We show how the probability distributions of a squared 
Bessel bridge and a squared Bessel process with or without absorption at zero are reduced to 
randomized gamma distributions. Moreover, for absorbing stochastic processes, we develop a 
new bridge sampling technique based on conditioning on the first hitting time at zero. Such an 
approach allows us to simplify simulation schemes. New methods are illustrated with pricing 
path-dependent options. 
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quasi-Monte Carlo method. 
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1. Introduction 

In this paper we study the exact path simulation of solvable continuous-time stochastic 
processes with transition probability density functions being obtainable in analytically 
closed-form. Despite the popularity of various approximation schemes for stochastic 
differential equations (SDEs), the precise path sampling of continuous-time Markov 
processes has certain advantages. Sampling from the exact probability distribution al- 
lows us to avoid introducing a bias and also to integrate along a path over an arbitrarily 
long time horizon. 

Our main motivation is the Monte Carlo pricing of path-dependent financial deriva- 
tives. The no-arbitrage price of a European-style option takes the form of a multi- 
dimensional integral along a path of an underlying asset price process. The usual 
procedure to the evaluation of such an integral is to employ the Monte Carlo method. 
Pricing of an American-style option reduces to solving a dynamic-programming prob- 
lem. Therefore, to apply the Monte Carlo method we have to sample paths from the 
exact distribution of the asset price process (e.g., see ifTTI ). 

More specifically, we study continuous-time Markov processes that arise from a 
squared Bessel (SQB) diffusion such as the squared radial Omstein-Uhlenbeck pro- 
cess (known also as the Cox-Ross-Ingersoll model), the constant-elasticity of diffusion 
model (with a power volatility function), and so-called hypergeometric diffusions ob- 
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tained from the squared Bessel process by means of a special combination of a change 
of measure and changes of variables (see ||4l|5l|6l). All these stochastic processes are 
broadly used in mathematical finance. Although for these models many fundamental 
quantities such as probability distributions of the first-hitting time at a barrier, maxi- 
mum and minimum values, and pricing formulas for barrier and lookback options can 
be obtained in closed-form, the Monte-Carlo method remains an important tool for the 
verification of analytical formulas and also for pricing Asian and American derivatives. 

As is shown in |;12|, the transition probability distributions of a squared Bessel pro- 
cess (without absorption at zero) and a squared Bessel bridge relate to the so-called 
randomized gamma distributions, which are mixture gamma distributions with a ran- 
dom rate parameter. The simulation of an SQB process with absorption at the origin 
is less studied in the literature. As is shown in H, the normalized transition density 
function of the SQB process is a gamma density which is randomized by a discrete 
probability distribution generated by a power series expansion of the lower incomplete 
gamma function. Therefore, to sample an increment of the random process we first 
simulate the absorption event and then sample from the normalized density function in 
case of surviving. Since we are able to derive the first-hitting time distribution of the 
SQB process with absorption at zero, it is possible to implement a completely different 
approach. First, we sample the first-hitting time, tq, at the origin. After that, we sam- 
ple the Bessel bridge with its value at time tq tied at zero. We show that the simplest 
realization of such an approach allows us to sample a path of the SQB process by only 
employing the gamma and Poisson probability distributions. 

The paper is organized as follows. Section 2 gives some basis results about the 
squared Bessel process and the squared Bessel bridge. Section 3 provides different 
sampling algorithms. In Section 4, we introduce other diffusion processes arising from 
the SQB process and provide simulation algorithms for them. Section 5 contains some 
numerical results. 



2, The Squared Bessel Process and Bessel Bridge 
2.1. The Squared Bessel Process 

Let us consider a Ao-dimensional squared Bessel (SQB) process {Xt)t>o obeying the 
stochastic differential equation (SDE) 

dXt = Xodt + vsfXtdWu Xtel= (0, oo), (2.1) 

with constant parameters Aq and > 0. The scale and speed densities are respectively 
s(x) = x~^~^ and m{x) = ^x'^, where ^ = ^ — 1 is called the index of the process. 
The left-hand boundary / = is entrance if ^ > 0, regular if — 1 < ^ < 0, or exit if 
H < —I. The right-hand boundary r = oo is natural. For the regular diffusion on I 
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the transition probability density function (PDF) is given by 

where /l = yu if / = is entrance or a regular reflecting boundary, and fl = \^\ if I = 
is exit or a regular killing boundary. 

For simplicity of presentation, we assume here that ly = 2. A simple scale transfor- 
mation xl""'^'"^ = Xl''°'^° \ X'q = (^) , allows us to modify without 



changing fi (i.e. ^' = ji"). 



2.2. The First Hitting Time Distribution 

In the case when / = is an absorbing boundary {fi < 0, p, = |/i|), the density in (12.21 ) 
does not satisfy probability conservation on J. The first hitting time (FHT), tq, at zero 
for the SQB process {Xt) starting at xo is defined by tq = inf{t : Xt = \ Xq = 
xq}. The PDF q{xo', t) for the FHT distribution is given by 

d 

q{xo;T) = -— p{T;xo,x)dx. (2.3) 
Jo 

By using that the transition PDF p satisfies Kolmogorov equations, we simplify the 
expression in j2.3l l to obtain 



1 d f p{t;xo,x)^ ^ 



5{x) dx \ m{x) 
As a result, we derive a closed-form expression for the FHT PDF: 



(2.4) 



A simple change of variable reduces the PDF in (I2.5l l to that of the gamma distribution 
G(a, /3) with shape parameter a = \p,\ and rate parameter fj = I. Therefore, the FHT, 
TO, can be sampled by using the formula tq = where Y ^ G{\fi\,l). 



2.3. The Squared Bessel Bridge 

Let < ii < t < t2- Consider a stochastic bridge generated by a continuous-time 
Markov process {Xt)t>Q G X with Xt^ and Xj, tied at x\ and X2, respectively. The 
bridge PDF 6 defined by h{t\,t2,t;x\,X2,x)dx = ¥{Xt E dx\xti = xi^Xt^ = X2} 
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can be expressed in terms of the transition PDF p of [Xi) as follows: 

. , X p{t-h;xi,x)p{t2-t;x,X2) 
b[ti , t2, t; xi,X2,x) = r . (2.6) 

p(t2 - tuXi,X2) 

Clearly, the bridge PDF h in ( |2.6l l integrates to unity thanks to the Chapman-Kolmogorov 
equation p{t2 — ti;xi,X2) = JjP{t — ti;xi,x)p{t2 — t; x, X2)dx. Notice that for the 
bridge density of a Gaussian process may also be derived in closed form by using a 
conditional multivariate normal distribution. 

The PDF of the squared Bessel bridge {Xt)Q<t<T conditional oxiXq = x and Xt = 
z is given by 

b{0, T, t; x,z,y) = ^, J , e"^ " f hiV^mi^JjT - t)) ^^.7) 



2t{T - 1) IfiiV^/T) 

where x = ^ ^ , y = and z = j^^^^ < t <T. 

Suppose that Xt is sampled conditionally on the FHT, T = tq. If t > tq, then set 
Xt = 0. Otherwise, if t < tq, we use the Bessel bridge with Xq and XT=Ta tied at x 
and z = 0, respectively. In the limiting case as z ^ 0+ in ( |2.7| |. we obtain 



5(0,r,.;.,0,,) = ^^^^|J exp(^-^J/,(V^). (2.8) 
Notice that the PDF in (I2.8ll has the same form as that in (12.2 



3, Simulation Algorithms 

In this section we present several algorithms for the precise path generation of the 
SQB process iXt). That is, for every time partition = to < ii < ■ ■ ■ < ^N, X > 
1, we sample a path-skeleton X = {Xq, Xi, . . . , X^), Xn = Xt„, from the exact 
multivariate probability distribution. The algorithms proposed below are all based on 
sampling from a randomized gamma distribution of the form G(a + Y, /?), where a + 
F > and /3 > are scale and rate parameters, respectively, and F is a nonnegative 
integer- valued random variable. As is mentioned above, we assume that ^ = 2, so all 
algorithms presented below deal with this case. In the general situation when v ^ 2, 
we proceed as follows. For given Aq, v, Xq, sample a path of the SQB process with 
/i = 2Ao jv^—X that starts at ( ^) ^ Xq by using one of algorithms in Figures[TJHl After 
that, rescale the path obtained by multiplying its values by (f )^. 



Exact Simulation of Bessel Diffusions 



5 



3.1. Randomized Gamma Distributions 



Suppose that a discrete random variable Y has discrete probabilities F{Y = n} = Pn, 
n = 0,1,2, ... . The PDF / of the mixture probability distribution G{a + Y, (3) admits 
the form of a series expansion: f{x) = Pn Y{a+n) ^ e~^-^ . 

Let us consider three choices for the randomizer Y of the gamma distribution G(a + 
Y, (3). The resulting distributions are called the randomized gamma distribution of the 
first, second, and third types, respectively. 

Let Yi P(A) be a Poisson random variable with mean A > 0. The randomized 
gamma distribution of the^zr^f type i?,G{Y\ + 6 + \, (3), 6 > - 1 , /? > 0, with the PDF 

/i(2/)=/3(^y y'/^e-^-^yie{Vmy), y > 0. (3.1) 

A discrete random variable Yz is said to have a Bessel probability distribution 
Bes{6, h) with parameters 9 > — \ and 6 > if 

r{Y2 = n}= \' ^.^.y n = 0,1,2,.... (3.2) 

Ie{h) n! Y{n + 9 + 1) 

This distribution is related to many other distributions, where the Bessel function / is 
involved in the density, including the squared Bessel bridge distribution (see lfT2l for 
details). The randomized gamma distribution of the second type is a mixture distri- 
bution G{Yi + 2Y2 + 9 + 1,(3), 13 > 0,9 > -I, where Fj ~ P((a + h)/{A(3)) and 
Y2 ^ 'Qes{9,^/ab/{2i3)) are independent Poisson and Bessel variates, respectively. 
For any positive numbers /3, a, b, and 9 > —1, the PDF is 

f2{y) = ^ ^ / e-^'^+'^"^-Pyie{V^)Ie{Vhl), V > 0. (3.3) 
l0{Vab/(2(3)) 

A discrete random variate Y^ is said to follow an incomplete Gamma probability 
distribution, which we simply denote by ir(^, A) with parameters A > and ^ > 0, if 

Notice that if ^ = 0, 1,2, . . ., then the distribution of is a truncated and shifted 
Poisson distribution thanks to the property 

7(m, a) , / x'^~^ 



. .1+X + ... + - — e ^, m = 0,l,2,., 

r(mj \ (m — Ij! 



We call a mixture Gamma distribution G{Y-i + \, (3),Y-i ^ IT{9, A), the randomized 
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input Xq > 0,0 = to < ti < • • • < tjv, M > -1 
for n from 1 to do 

Xn-l 



1 



^ "^(j^n ^n— l) 

x„~g(^f„ + ^ + i, 

end for 

return (Xo,Xi, . . .,Xn) 



in—i) 



Figure 1. The sequential sampling method for modeling an SQB process without absorption. 



gamma distribution of the third type. The PDF is 



My)=P 



T{e) (p 



7(0,A) VA 



-e/2 



(3.5) 



3.2. Simulation of Processes witliout Absorption 

The randomized distribution of the first type is closely connected with the transition 
distribution of a squared Bessel process [Xt) without absorption (i.e. > 0, or /i G 
( — 1,0) and a; = is a reflecting boundary). The conditional distribution of Xt, t > 0, 
given Xq = xo > 0, is then a randomized gamma distribution of the first type. The 
transition PDF in (I2.2| | with v = 2 has the form of the PDF f\ in (13.11 1 with 9 = p, 
f3 = I /2t, and A = xo/2t. Therefore, we have the following sampling scheme: 



Xtr^G{^i + Y + l,l/2t), where y - P(xo/2t), t > 0. 



(3.6) 



The sampling algorithm is presented in Figure [T] 

A path of the standard squared Bessel bridge can be generated using the second type 
randomized gamma distribution. The bridge PDF in i2Jl reduces to that in ( 13.31) by 
setting a = x/t^, h = z/{T — t)^, j3 = 2t{T-t) ' ^^'^ ^ ~ Th^ti' -^t conditional 
on Xq = X and Xt = z, < t < T, can be obtained by generating two independent 



random variables Y 



Xt^G{Y + 2Z + ^i + l 



j_ 

2T 



and Z ^ Bes ( fi, 



' 2t(T-t) 



)• 



and then 



3.3. Sequential Simulation of Processes with Absorption 

Assume that a stochastic process {Xt)t>Q G R+ admits absorption at the origin. 
For example, for an SQB process we have that ^ < and a: = is a killing boundary 
or exit. Clearly, the transition PDF p given by ( |2.2| | with fi = /i < 0, does not 
integrate to one. Let us define the probability Pg of surviving before time t and the 
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input Xq > 0,0 = to < ti < • • • < ^Af, M < 

To ^ OO 

for n from 1 to do 
if To = oo tlien 

Un^V (0,1) 

if Un < Pa tlien To ^ tn 

end if 

if t„ < To tlien 

^■~'K''''' 2(J-tl-,) ) 

else 

Xn^O 

end if 
end for 

return (Xo, Xi, . . . , X^) and fo 
Figure 2. Tlie sequential sampling method for an SQB process with absorption at the origin. 



probability Pa of absorption before time t for the process (Xt) started at Xq = x: 

Ps{x; t) = p{t; X, y)dy > and Pa{x; t) = I - Ps{x; t) > 0. 
Jo 

Observe that the actual transition probability distribution is then a mixture of continu- 
ous and discrete probability distributions with the following generalized PDF: 

p{Xo - X,) = P,{Xo;t) ■ ( ^p^^^^.y ) + Pa{Xo;t) ■ 5{X,), 

where 5 denotes a delta function. 

By using (12.51) . we obtain the following probabilities of surviving and absorption of 
the SQB process before time t: 

P.{x;t) = P{to >t} = iMil and P,(x;i) = P{to <t} = 

where 7(0, x) and r(a, x) are the lower and upper incomplete gamma functions, re- 
spectively. The normalized transition PDF of the SQB process conditioned on the 
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input Xq > 0, = to < ti < ■ ■ ■ < In , fj, < 

F~G(H,1), To^^ 
for n from 1 to do 



if tn < To tlien 



Xn-l{TQ — tn) 



2(to — tn-l){tn — tn-l) , 



1, 



(to — tn){tn — tn-l) 

else 

end if 
end for 

return {Xq, Xi, . . . , X^) and tq 



Figure 3. The sequential sampling method conditional on the FHT, tq, for modeling an SQB 
process with absorption at the origin. 



survival of the process before time t is 

p{t-x,y) _ r(lMl) {x\'^ e-("+^°)/2* {V^o\ 

As is seen, the function in the right-hand side of (13.71) reduces to the form of ( 13.51 1 
with 6 = \^\, A = x/2t, and fj = l/2t. Thus, the above normalized transition PDF 
follows the randomized gamma distribution of the third kind G(F + 1, l/2t), where 
Y ~ ir{\^\,x/2t). As a result, we obtain the sampling algorithm given in Figure[2l 
The algorithm returns a sample path X and an approximation, tq e [ti, . . . , ^at, oo}, 
of the FHT, tq. 



3.4. Bridge Simulation of Processes with Absorption 

Consider again the SQB process (Xt) with absorption at the origin. Since the first 
hitting time PDF q{xQ; r) is available, we may first sample the FHT, tq, and then 
simulate a path of {Xt)t>o conditional on tq by using the bridge distribution. As is 
seen from (12. St . the PDF of Xj, < t < tq, conditional on Xq = x and Xr^ = 
is reduced to the PDF /i in (13. Il l of the randomized gamma distribution of the first 
type with 6 = A = '^^J""*^ , and {3 = jtirl-t) • ^ result, we obtain a sequential 
sampling algorithm conditional on the FHT (see Figure[3l). 

At last, in Figure |4l we provide the full bridge sampling algorithm, where a path 
X = (Xq, Xi , . . . , Xtv), N = 2'^, k > 1, is sampled at the time points in the following 
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order of generation: 

*Af)^Af/2)^JV/4)^3Af/4>iAr/8,i3A'/8i^57V/8,^7Ar/8, • • • , ^2, ^6, • • • , ^1 , ^3) ■ ■ ■ ,^JV-1 



Here, we use that the bridge PDF in (I2.7l l with /i = \fj,\ reduces to that in (13.31 1 by 
setting a = x/t^, b = z/{T — t)-^, f3 = 2t(T-t) ' ^ ~ 1^1- Such abridge samphng 
algorithm is very useful for the quasi-Monte Carlo pricing of path-dependent options. 



4, Generating Paths of the CIR, CEV, and Hypergeometric Diffusions 
4.1. The CIR Process 

Consider the Cox-Ingerssol-Ross (CIR) diffusion process {Yi)t>Q G X = i?+ solving 
the SDE 

dYt = {\o-MYt)dt + v^tdWt, (4.1) 

where constant parameters Aq, Ai, and v > 0. The respective scale and speed den- 
sities are s(x) = x~^^~^e'^^ and m(x) = ■^x'^e~'^^, where k = The boundary 
classification of the CIR process is equivalent that of the SQB process. For the regular 
diffusion on X, the transition PDF is 

p{t;x,y) = cte^^' i^j e'^^^-^ '*+-)/^ (2c^^^^j , (4.2) 



where Cf = K/{e^^* — 1) and jl is defined as for the SQB process in Section|2l 

The CIR process is reduced to an SQB process with the same parameters Ao and f 



by means of scale and time transformation, Yf = e ^^^Xg u\, where the monotonic 



time-transformation function is defined by 

it if Ai = 0, 

I ^ io,^o. 

The transition PDF for the CIR process relates to that of the SQB process as follows: 

p^^^^\t;x,y) = e^''p'^^^^\sx,{t);x,e^''y). 

If a reflecting boundary condition is imposed at a; = 0, or the origin is entrance, 
then the CIR diffusion is a conservative stochastic process. The corresponding tran- 
sition density is given by (I4.2| | with jl = ^ > —I. The transition distribution of 
the conservative CIR model reduces to the randomized gamma distribution of the first 
type. The respective SQB process admits no absorption at zero and can be simulated 
by the sequential method in Figure [T] 

Consider the case where x = is a killing boundary or exit, so the transition PDF 
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input XQ>0,0 = to<ti <■■■ <tN,N = 2'',^i<0 
F~G(|/.|,1), ro^^ 
if ijv <To then 



2ToiAr y' \ 'i7v(To-iAf) 

else 
end if 

for I from 1 to /c do 

for 771 from 1 to 2'"^ do 
71 = (2m - 1)2^-' 
if tn > To then 

else 

if tnj > then 

Xrn {to — tn) 



Yr 



2(to — tnf){tn — tni) 

X„ ~ G (V„ + ImI + 1 ~ 



"1 



' {to — tn){tn — tni) , 

else 

Xni{tn2 tn) ^712 (^^ ^^1) 



Bes ( |Ai| 



2(to — tnj){tn — irii, 



X„~G(y„ + 2Z„ + |M| + l, 



end if 
end if 
end for 
end for 

return (Xq, Xi , . . . , X^) and tq 



2(tyi tm){tfi2 ^n) 



Figure 4. The full bridge sampling method conditional on the FHT, to, for modeling an SQB 
process with absorption at the origin. 
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is given by (14 .21 1 with p. = \^\, where /i < 0. The FHT, tq, at zero for the CIR model 
is given by 

^ inf{t : = 0} ^ inf{i : X,,^(,) = 0} ^ ^xUrr^'^), 

where we define s^^{t) = cxd if t > sa, (oo). The corresponding PDF is given by 

q^^^^\xo;T) = e^i^(7(^Q^)(xo;sA,(T)). We have that P{r^^^^' < oo} = P{Tf < 
SA,(oo)}. 

Clearly, the sampling of a CIR path at times ti, i = 0,1, ... ,N, reduces to the 
sampling of an SQB trajectory. The method for sampling a path and the FHT, tq, is 
given as follows. 

Step 1. Set times Si = s\^ (ti), i = 0,1, . . . , N. 

Step 2. Obtain a sample path {Xq, Xi, . . . , X^) of the SQB process at times Si, i = 
0,1, . . . , N, and the FHT, Tq^'^^\ (or its approximation fo) by using one of 
the algorithms in Figures [TJHl 

Step 3. Set Yi = e'^'^'Xi, i = 0,l,...,N. 

( _w {SQB). .f {SQB) ^ , s 
{CIR) _ I S (Tq ) if < SAi(oo) 



Step 4. SetTQ 

Step 5. Return (Fq, Yi , . . . , Fat) and 



oo otherwise 

(CIR) 



4.2. The CEV Diffusion Model 

The constant elasticity of variance (CEV) diffusion process {Ft}t>o obeys the stocha- 
stic differential equation dFt = rFtdt + 5F^^^ dWt, t > 0, Fq > 0, where r, 5, f3 are 
real parameters. We assume here that S > and /3 < 0. 

The boundary F = of the state space (0, oo) is regular if (3 < —0.5 or exit if 
—0.5 < P < 0. Here we consider the case where the endpoint F = is a killing 
boundary. The transition PDF po{t; Fq, F), Fq, F > 0, t > 0, for the CEV process 
(Fj"') with zero drift (r = 0) takes the form 

The density po{t; Fq, F) does not integrate (with respect to F) to unity for t > 0, since 
F = is an absorbing point. 

Ir) (0) 

A drifted CEV process F^ with r 7^ is obtained from F^ by means of scale 
and time transformation: f/^' = e^^F^^K.., where sx, is given by (I4.3l l with Ai = 
2rj3. The resulting transition density p^. with r 7^ is given by pr{t; Fq, F) = 
e-'*po{sx,{ty,FQ,e--'F). 
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The Monte Carlo simulation of the CEV diffusion is based on the reduction of it 
to the CIR or SQB process by using the mapping X(F) = ^r^- There are two dual 
approaches: 

(i) First, eliminate the drift and then, by using the mapping X, reduce the driftless 
CEV process to an SQB process defined by Xt = X(Fj°'), t > 0, with Aq = 
2+1/(3 and v = 2. Sample a path of the SQB process and then obtain a path of 
the driftless CEV process by applying the mapping F(a;) = {S^(3^x)~^^^^. After 
that, restore the drift using the time and scale transformation. 

(ii) By using the mapping X, reduce the drifted CEV process to a CIR process defined 
by Yt = X(f/''^), with Aq = 2 + 1//3, Ai = 2r/3, and u = 2. The resulting 
CIR process can be obtained from an SQB process by means of time and scale 
transformation. Sample a path of the CIR process and then obtain a path of the 
CEV model by applying the inverse mapping F. 

The FHT, tq, at zero for the CEV diffusion model is given by 

rf . inf{. : = 0} i rf I .."/....(rf 

Notice that if a reflecting boundary condition is imposed at F = when (3 < —0.5 
(or /3 > and hence F = is entrance), then the CEV diffusion is a conservative 
stochastic process. The corresponding transition density (for the case with (3 < —0.5) 
is given by (I4.4| | with the replacement Ij_ . By analogy with the CIR model 

2\0\ 20 

without absorption at zero, the transition distribution of the conservative CEV model 
reduces to the randomized gamma distribution of the first type, hence the algorithm in 
Figure [T]is applied. 

4.3. Diffusion Canonical Transformation 

Several families of analytically solvable diffusions can be derived from known under- 
lying diffusion processes. We refer to this construction as the "diffusion canonical 
transformation" methodology (see ||4j|5]|6l for details). 

Let us start with a one-dimensional time-homogeneous regular diffusion {Xt)t>o G 
X = (/,r), — oo < Z < r < oo, defined by its infinitesimal generator: {Q f){x) = 
^v^{x)f"{x) + X{x)f'{x). The functions A and ly denote, respectively, the (infinitesi- 
mal) drift and diffusion coefficients of the process. Consider two linearly independent 
fundamental solutions ip'^ and ip~ of the differential equation {Q ip){x) = sip{x), 
s e C, X e X, such that for real values s = p > the solutions and (p~ are 
respectively increasing and decreasing functions of x (see, e.g., jSj). 

Let us introduce another diffusion (X^ ''')t>o G X with generator 

1 / u' (x)\ 

(a^"' /)(x) ^ -i^\x)f{x) + Mx) + ,^\x)^ fix) , (4.5) 

2 V Up{x)J 
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where a strictly positive function Up{x), p > 0, is a linear combination of Lpp. Up[x) = 

q\ip'^{x) + q2(p~{x), (71,2 > 0, gi + 92 > 0. A transition density p^p for the X^^^- 
diffusion is then related to a transition density px for the X-diffusion as follows: 

pP{t;xo,x) = e-P'^^^^Px{t;xo,x), x,xo&T,t>0. (4.6) 

Now we consider an F-diffusion [Ft = F(X('''), t > 0} defined by strictly mono- 
tonic real-valued mapping F = F(x) with F', F" continuous on X and having in- 
finitesimal generator {gfh){F) = \a^{F)h"{F) + rFh'{F), where F e Jp = 
(min{F(/+), F(r— )}, max{F(/+), F(r— )}), and r is a real constant so that p + r > 0. 

The transition PDF pp for an F-diffusion {Ft)t>o is related to the transition PDF 
for the underlying X (or X'^p^) diffusion as follows: 

a{F) Mp(X(Fo)) 
Here X = F~Ms the inverse map. F admits the general quotient form: 

C ^ _ Cl</?p+r(2;) + C2'/'p+r(2;) _ Vp+r{x) 



qiip+{x)+q2ipp{x) Up{x) 



(4.8) 



where ci and C2 are real constants. For a full classification of strictly monotonic maps 
of the form (14.81) see JU. The diffusion coefficient function is 

^(^)^K:^)|H:M^ ^^X(F), F^Tp, (4.9) 



u 



(x) 



where we define the Wronskian W{x) = Up{x)v'p^j.{x) — u'p{x)vp+r{x) ■ 

In the next two subsections we present two examples of hypergeometric diffusions. 
The concluding subsections gives a general simulation algorithm. 



4.4. The Bessel- K Diffusions 

Here we specifically consider a 4-parameter Bessel K-family arising from an under- 
lying (Ao-dimensional) squared Bessel process with a positive index ^. We use the 
generating function Up{x) = tpp{x) = x~^^^Kp^ {2^/2px/v) and the mapping: 

{2^2{p + r)x/v) 

F{x) = c — ^ — — ^ , , ^ , (4.10) 

Kp [l^lpxlv] 
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where c, p, u, and fi are independently adjustable positive parameters, and r > —p is a 
real constant. The functions / and K denote the modified Bessel functions of the first 
and second kind, respectively, (see [Tl for definitions and properties). 

The function F(a;) (and the respective inverse X(F)) maps x e (0, cxd) and F e 
(0, cxd) into one another. The transformation ( I4.10l l hence leads to a family of processes 
(Ft) e (0, oo) with the diffusion coefficient function 

a[t(x))-cvzy Kl{l^) + K,{1^) J ^ ' 

Lemma 4.1 (Campolieti and Makarov, fU^j). The processes of the Bessel K-family 
obeying the SDE dFt = rFtdt + CF{Ft)dWt with ( l?JOl )- d?T/D have the following 
boundary classification: the boundary F = is exit if fi > I or is a regular killing 
boundary ifO < n < 1; the boundary F = oo is non-attracting natural. Moreover, 
the discounted process (e~'^*i^t)t>o is a martingale. The transition PDF pp is given 
by 114. 7i with ^{x) = Vy/x, and a and px respectively specified by A4.1U and A2.2h 

The density, q{Fo;T), for the FHT at the origin for a Bessel-K process started at 
Fo > is readily derived by using equation ( 12.4b . giving the generalized inverse Gaus- 
sian distribution: 

q{Fo-T) = J^^^^^—^r-^-^e-i^-^--l^'\ T > 0, xo = X(Fo). (4.12) 
4.5. The Confluent-ZY Diffusions 

The confluent hypergeometric family of F-diffusions arises from an underlying CIR 
process with p. > 0. Here we specialize to the confluent-U family with generating 
function Up{x) = ip~ {x) =U{v, fx+ 1,kx) and mapping 

M{v + ^,fi+ 1,KX) 

F(x) = c ^ (4.13) 

U[v,p. + 1, KX) 

where v = p. = ^ — 1, k = and c are arbitrary positive constants, and r > 
—p. The confluent hypergeometric functions Ai and U are two linearly independent 
solutions to Kummer's differential equation (see HI for definitions and properties). 

The function F(z) maps x e (0, oo) onto F S (0, oo) and is monotonically in- 
creasing. This transformation leads to a family of processes {Ft) e (0, oo) with the 
diffusion coefficient function 

a{f{x)) = CKV^ u^-[^,i.U,.x) + p.l4) 
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Lemma 4.2 (Campolieti and Makarov, ||4l|6l). The processes of the confluent U -family 
solving the SDE dFt = rFtdt + a{Ft)dWt with i4l3l - i[4l4\l have the same boundary 
classiflcation as that for the Bessel-K in Lemma \4J] Moreover, the discounted process 
{e~^* Ft)t>Q is a martingale. The transition PDF pp is given by l\4. 71 1 with i'{x) = 
v^fx, and a and px respectively specifled by 0. 14^ and f |4.2D . 

The density for the first-hitting time at the origin, q{FQ\ r), for a confluent-i/ process 
started at Fq > is 

-KxaT{T)fn-l\\v-\fi,n-f \Y-v 

q[F,;r) = \T\r)\^- J ^ i\ \r/ ' " > ^^.IS) 

' ' U{v,ii+1,k.xq)Y{v) 

where xq = X(Fo) and we use the time change T(r) = -x^r ' "^^^ latter func- 
tion in J4.15I I is known as a Tricomi exponential PDF (see 111) given by p{T) = 

^-zTq-a-lM _|_'2-\fe-a-l 

-r-r-—, — ; — ^ , T > 0, where a = v,h = + \, z = kxq,. It integrates to 

T[a)U(a, 0, z) 

unity thanks to the integral representation of U (see HI). 
4.6. Simulation of F-Diffusions 

We generalize the sampling algorithms for an SQB process presented in Figure [3] and 
Figure |4l Within that approach a path is sampled conditionally on the FHT at zero. 
The Bessel-K and confluent-W diffusion models are both absorbing at zero and have 
the first-hitting time distribution in analytically closed-form. For a sampling algorithm 
we only need to obtain the distribution of the respective bridge process. In doing, 
so we use one important observation that the distribution of an F-diffusion bridge is 
reduced to the distribution of a bridge of the respective underlying diffusion (e.g. the 
Bessel and CIR bridges). 

By applying the analogue of formula (I2.7| | for an F-diffusion with PDF pF{t; Fq, F) 
in place of the PDF p{t; x, y), and using the representation (14.7b . we have the following 
expression for the bridge PDF of an F-diffusion with Ft^ and Fj^ tied at Fi and F2 
respectively: 

bF{tut2,t;FuF2,F) = ^l2i^5|^)(ii,i2,i;X(Fi),X(F2),X(F)) (4.16) 
= '^^^^bx{h,t2,t;X{F,),X{F2),X{F)) 

where hx and h^^ denote the bridge PDFs of the diffusions [Xt) and {x[''''), respec- 
tively. Here, after plugging ( I4.7| | in the formula of the bridge PDF 6/^, we first cancel 
Jacobians - and then cancel Doob's factors of the form e~^* "''frj . If follows from 

cr Up(x) 

( 14.161 ) that an F-diffusion bridge is obtained by applying the mapping function F to 
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the bridge process for the underlying diffusion {Xt) with Xf, and Xt^ tied at X(Fi) 
and X(F2) respectively. For example, in the particular case of the Bessel-K diffusion 
when the underlying process {Xt) is a squared Bessel process, the F-bridge is just a 
nonlinear transformation of a standard Bessel bridge. 

Our primary goal is to sample a path skeleton (Fq, Fi, . . . , F/v), Fi = Ft- of an 
F-diffusion at times ti, i = 0,1, ... ,N, = to < ti <■■■< In = T, for a given 
initial condition Ft=o = Fq. The simulation scheme based on the bridge distribution 
is as follows: 

Step 1. Sample the FHT, tq, from the GIG or exponential Tricomi distribution for the 
Bessel-K or Confiuent-W model, respectively. 



Step 2. Obtain a sample path (Xq, Xi, . . . , Xj^) of the respective underlying process 
(the SQB or CIR diffusion) conditional on Xq = X(Fo) and Xr^ = 0. 



Step 3. Apply the respective mapping function F to obtain a sample path of the F- 
diffusion model: Fi = F(Xi), i = \,2, . . . ,N . 

The main result is that this simulation scheme allows us to avoid a direct sampling 
from complicated transition probability distributions. 

Let us present an alternative approach from ||5l to computing mathematical expec- 
tations of path functionals of the form Q = E[/(Fi, F2, . . . , F/v)|Fo] for F-diffusion. 
By using a path integral approach, the expected value of such a path functional can be 
represented as a multivariate integral: 



Xk)dxi ■ ■ ■ dxN- 



Q = [ f{f{xi),...HxN))e-''^^^^^l\px{tk-tk-uXk-u 

Jrn Up{xo) ^^Jj 

The integral above may be estimated by the Monte Carlo method. The underlying 
diffusion is simulated by sampling from the exact transition probabilty distribution. 
The resulting unbiased estimator of the path integral Q takes the form: 

C = /(F(XO,...F(X^))e-^^^^^, 

Up[XQ) 

where the path {Xo,Xi, . . . ,Xjv) is sampled by using one of the algorithms from 
Section [3] Notice that we cannot use the Euler method (or any other approximation 
method, which does not guarantee the positiveness of the approximation process) since 
the estimator £^ is infinite if Xyy = 0. Using large and small argument asymptotics of 
the Bessel function K and Kummer function U, we obtain that the variance of ^ is 
finite if /i < 1. Notice that the use of an exact simulation method allows us to lift this 
restriction. 
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5, Simulation Study 

5,1, Simulation of Randomizers 

The three discrete probability distributions used in the construction of randomized 
gamma distributions are all log-concave and unimodal as is stated below. 

Lemma 5,1, Let Y be a Poisson, Bessel or incomplete gamma random variable. The 
distribution ofY is log-concave. That is, the ratio ¥{Y = n + 1}/¥{Y = n} is 
decreasing in n. Furthermore, the distribution of Y is unimodal and has a unique 
mode or two modes at consecutive integers. Moreover, one mode is always located 
at m = [AJ for the Poisson distribution, at m = + 6*2 -0)/2\ for the Bessel 

distribution and at m = max(0, [A — ) for the incomplete gamma distribution. 

Proof. See [S] ID for the proof for the Bessel and incomplete gamma distributions, 
respectively. □ 

To generate a Bessel or incomplete gamma random variate, we can use a generic 
acceptance-rejection (A-R) method from IT) stated below without proof. 

Lemma 5,2 (Devroye, Q). For any discrete log-concave distributions with mode at 
m, we have, for all n > 0: Pn < Pm min 1 1, gi-Pm\n-m\ | _ 

As an alternative sampling method we use the inversion method by chop-down 
search (C-D-S) from the mode m. Such a sampling method for a discrete distribu- 
tion with probabilities {pk}k>o is based on the numerical inversion of the CDF F by 
the formula F~\u) = argminjn > | u — '^^^aPk < 0}, u e [0,1]. 

It is well known that the computational cost of such a method has the lowest possible 
value if and only if the vector of discrete probabilities is arranged in increasing order. 
Instead of the preliminary computation of probabilities followed by sorting of them, we 
start the search algorithm at the mode m and then successively calculate probabilities 
of values to the left and to the right of the mode choosing the largest one. Notice 
that probability pm need only be computed once, and that other probabilities can be 
obtained by using simple recurrences. 

We now present some numerical results comparing the two methods of simulation 
of P(A), Bes(6',b), and ir(6', A) random variables. For each of the two methods, one 
million values are sampled. For simulation of each of the Poisson random variables, 
the parameter A is allowed to vary as a continuous uniform random variable. For 
the two parameter Bessel and incomplete gamma distributions, the first parameter is 
allowed to vary as a continuous uniform random variable while the second parameter 
is held constant. Then the procedure is repeated by allowing the second parameter to 
vary while the first one is held constant. Results of these tests are given in Tabled] 

Table 1 shows that the chop-down search method from the mode is significantly 
faster than the acceptance-rejection technique for generating random variables in ev- 
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Table 1. Comparison of the acceptance-rejection and chop-down search methods for 
the Poisson, Bessel, and incomplete gamma distributions. 

Distribution A-R Method C-D-S Method 







Time 


No. of Iter. 


Time 


No. of Iter 


P(A) 


A ~ U(0, 1000) 


189.9 


2.6 


35.2 


34.2 


Bes(6',6) 


e ~ u(o, 1000), b = 10 


220.6 


1.6 


100.4 


1.1 


Bes{6l,b) 


9 = 10, 6 ~ U{0, 1000) 


414.1 


4.0 


103.3 


17.3 


ir(6', A) 


e ~ u(o, 100), A =10 


336.6 


3.7 


51.1 


1.8 


ir{e,\) 


6» = 10,A ~ U(0, 1000) 


363.7 


3.9 


51.4 


10.8 



Note: Time in seconds and average number of iterations for the simulation of 10^ random variables 
from the Poisson, Bessel, and incomplete gamma distributions using the acceptance-rejection (A-R) and 
chop-down search (C-D-S) methods. 

ery case and is a much better choice for simulating random variables when it can be 
implemented. 

5.2. Comparison of Sampling Schemes for the SQB Process 

In this section we aim to compare the following three sampling schemes. 

1) Sequential sampling conditional on the FHT tq with the use of the randomized 
gamma distribution of the first kind. 

2) Bridge sampling conditional on the FHT tq with the use of the randomized gamma 
distribution of the second kind. 

3) Unconditional sequential sampling with the use of the randomized gamma distri- 
bution of the third kind. 

We start by sampling multiple paths of the SQB process over a discretized partition of 
a time interval [0, T], = io < < ■ ■ ■ < tjy = T, using one of the three methods 
just mentioned. Then we average these sample paths in order to approximate the mean 
of the SQB process. To study the sampling algorithms, we compare our sample means 
to the true mean of the SQB process as well as the time required to simulate a set 
number of sample paths of the process. 

For calculation of the mean of the SQB process, we use the formula 

which is valid for ^ < and u = 2. The expression is derived by considering the mo- 
ment generating function of the SQB process at time t and using the small asymptotics 
of the Bessel / function. 

To this end, we look at the largest amount by which the sample mean, fit, dif- 
fers from the true mean, /it = E[Xi], (the maximum absolute error) at times ti,i = 
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0, 1, . . . , A^, given by max \fit^ — /ijj. We will also examine the largest sample 

1=0,1,..., A'' ' ' 

Standard deviation of the process, given by max d-tj\/n, where n is the sample 

size. After simulating one million sample paths for each of the three sample schemes 
and averaging them, we obtain the data shown in Table |2l From this data, we can see 

Table 2. Comparison of sampling schemes for the SQB process. 







Scheme 1 






Scheme 2 






Scheme 3 






Time 


MAE 


MST 


Time 


MAE 


MST 


Time 


MAE 


MST 


-0.25 


1741 


.00244 


.00271 


5044 


.00297 


.00270 


2501 


.00534 


.00270 


-0.5 


1600 


.00111 


.00252 


4462 


.00191 


.00251 


2280 


.00158 


.00251 


-1.5 


953. 


.00193 


.00144 


2614 


.00240 


.00145 


1406 


.00149 


.00145 



Note: Time in seconds, maximum absolute enor (MAE), and the maximum standard deviation (MST) 
taken from the average of 10* sample paths of the SQB process using sampling schemes 1, 2, and 3 
respectively for varying values of /i. For all three choices of /i, we set Xo = I, T = \, v = 2, and the 
partition of [0, T] to be 0, 3^ , ^ , . . . , 1 . 

that sampling scheme 1 is the fastest one. Scheme 2 is much slower than schemes 1 
and 3 since it involves sampling from the Bessel distribution. 

5.3. Sampling from the GIG and Tricomi Exponential Distributions 

A common approach to sampling from a nonstandard probability distribution is to use 
an acceptance-rejection method. This approach is employed in jj) and ||9l for sampling 
from the GIG and Tricomi exponential distributions, respectively. If the parameters of 
a probability distribution remain constant, then a much faster sampling technique is the 
one that is based on the numerical inversion of a distribution function. To sample from 
a continuous CDF F by using the inverse transform method, we generate a uniformly 
distributed on (0, 1) random variable U and then set X = F~\U), where F~' the 
inverse of F. 

In cases where the inverse of F can not be expressed in closed-form, the inverse 
transform relies on numerical approximation. A root-finding method such as Newton's 
method or the bisection method can be applied to solve equation F{X) = U, U & 
(0, 1). A faster approach is to compute the CDF on a fine mesh and then approximate 
the inverse of the CDF by some simpler functions. The simplest method is to use a 
piece-wise linear interpolation. In ifTOl a fast and efficient variate generation method 
is proposed. In that method, the inverse CDF F~^ is approximated by the Hermite 
interpolation functions. For a given partition I = xq < xi < ■ ■ ■ < Xn = r of the 
support (/, r) of a CDF F, the distribution function is computed by either integrating 
the density function on each subinterval {xk-i,Xk), or by employing an ODE solver, 
since the CDF F solves a simple ODE F'{x) = f{x), x e (^t'), F(l) = 0, where / 
is the respective PDF. 
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5.4. Path-Dependent Options 

This section reviews some discretely-monitored path-dependent options that will be 
used for pricing options in the following subsection. First, we assume that we have 
sampled a path of a asset price process (Ft) over a discrete time partition, T = 
{ij}j=o,i,...,JV> of the time interval [0,T], T > 0. Let the values of process {Ft) at 
time points t = Ube denoted by F^, for alH = 0, 1 , . . . , N. 

The payoff function of an Asian-style option depends on the arithmetic average of 
the underlying asset values: An = jj S^li Fi. For an average price call option, the 
payoff to the option holder at time T is {Am — K)+ where K is the strike price and 
(a;)+ = max(x,0). The average price put option is defined similarly. Its payoff at 
time r is {K - An)+- 

The second type of path-dependent options we will price are lookback options. In 
this case, the payoff functions depend on the maximum, = max^F^, or the 

minimum, rtiM = min Fi, values of the underlying asset price attained during the 

i=0,\,...,N 

option's life, [0, T]. A standard lookback call gives the right to buy at the lowest price 
recorded during the options life. Hence, the payoff to the holder at time T is FM — rtiM- 
A standard lookback put gives the right to sell at the highest price recorded during the 
options life. Thus, the payoff at time T is Mn — Fn- 

5.5. Pricing Patli-dependent Options under Nonlinear Volatility Models 

In this section we present some numerical results regarding pricing Asian and look- 
back options under the CEV, Bessel-K and Confluent-ZY families of diffusions using 
Monte-Carlo algorithms based on generating from randomized Gamma distributions. 
Specifically, we look at a plain sequential Monte-Carlo sampling method (MCM) and a 
randomized quasi Monte-Carlo method (RQMCM) which uses digital scrambling via 
a Sobol's sequence for the randomization. For the Bessel-K and Confluent-W models, 
we also use the weighted method (MCMW) described in Subsection l4.6l One million 
simulations are completed for each payoff function and are then averaged to get the 
final option pricing results. For the RQMC method, these 10^ simulations correspond 
to 100 randomizations and 10000 simulations per randomization. 

In the tests that follow, we fix the value of the annual local volatility function 
o'ioc('S'o) = 0-25 at the initial asset price = 100. The strike price is K = 100. 
The interest rate is r = 0.02 per annum and all options have six months to expiration: 
T = 0.5. The number of asset price observations is = 128. First we look at pricing 
under the CEV model. For the CEV model, (Jioc{Sq) = SS^. Typical observed values 
of the CEV elasticity parameter /3 are strongly negative so we choose /3 = —2. Then 
we choose the parameter 5 so that it satisfies SF^ = 0.25. This yields 5 = 2500. Next 
we consider the Bessel-K subfamily of diffusions. To ensure that aioc{Fo) = 0.25 the 
following parameters are used: p = 0.001, r = 0.02, c = 154.4870, // = 0.25, and 
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V = 2. The last pricing model considered here is the Confluent-W family of diffusions. 
Specifically, we examine the case where c = 788.3679, p = 0.001, Ai = 0.0009, 
p, = 0.25, and v = 2. Table [3] contains option pricing results corresponding to these 
models. The prices reported are obtained using the RQMC method. Table |4] reports 
the computational cost of pricing the average price Asian call using the three methods. 

Table 3. Pricing path-dependent options under the three models using the RQMC 
method. The value of the sample standard error is given after the ± sign. 

Model Asian Call Asian Put Lookback Call Lookback Put 

CEV 4.30237±.00081 3.80260±.00160 14.55220±.00255 12.09087±.00300 

Bessel-K 4.28605±.00049 3.797 17±.00033 13.15557±.001 13 13.23640±.00081 

Confluent-W 4.28724±.00049 3.79922±.00032 13.31 158±.00093 13.11594±.00084 



Table 4. Computational cost of pricing the average price Asian call option. 



Model 


Method 


Smpl.Var., a' 


Time (sec) 


Cost, a-T 


Relat. Cost 


CEV 


MCM 


32.574 


7438 


242296 


52.4 




RQMCM 


0.065 


70762 


4622 


1.0 


Bessel-K 


MCMW 


33.044 


33291 


1100088 


52.6 




MCM 


41.830 


10506 


439444 


21.0 




RQMCM 


0.235 


89029 


20895 


1.0 


Confluent-W 


MCMW 


31.636 


33312 


1053853 


49.4 




MCM 


40.801 


10174 


415122 


19.5 




RQMCM 


0.238 


89715 


21308 


1.0 



As seen in Table IH the RQMC method offers a clear improvement in reductive 
cost over the plain MC method. On the other hand, the weighted method offers no 
improvement in cost at all, mostly due to its relatively large computational time. The 
extra time required for the weighted method is partly due to the computation of special 
functions in the weight. It could also be attributed to sampling more points in each 
of the sample paths for a price process {Ft). When conditioning on the FHT tq and 
sampling at time t, we check first whether t > tq. \f t > tq we do not have to sample 
from any probability distributions since Ft = 0. When using the weighted method 
we are looking at the case with no absorption so we don't have this benefit. In other 
words, for every point of the discretized sample path, we must sample from probability 
distributions which takes up more time. This combined with the fact that fox p> 1 we 
have no guarantee that the mean of the weighted estimator is finite makes the weighted 
method a poor choice for pricing options. We have a much better choice in the exact 
sampling method. 
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